library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=3)

model.invalid <- "
data{
int N;
int K;

int D;
int<lower=1, upper=D> Dist[N];

vector[N] y;

matrix[N,K] X;
}

parameters{
real alpha;
vector[K] beta;
vector[D] nu_d;

real<lower=0, upper=10> sigma_y;
real<lower=0, upper=10> sigma_d;
}

model{
beta ~ normal(0, 0.1);

sigma_y ~ uniform(0, 10);
sigma_d ~ uniform(0, 10);

y ~ normal(alpha + X*beta + nu_d[Dist], sigma_y);

nu_d ~ normal(0, sigma_d);
}
"

stanModel.invalid <- stan_model(model_code=model.invalid)




load("JHRED2014_Invalid.RData")




###
### sntv invalid
###

# election-year matrix
year_list <- sort(unique(sntv_district$year))
year.mat.sntv <- matrix(0, nrow=nrow(sntv_district), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.sntv[sntv_district$year==year_list[i],i-1] <- 1
}

x.mat.sntv <- cbind(sntv_district$kakusu.mean, log(sntv_district$ku_ncand), 
                    sntv_district$ku_m, year.mat.sntv)

sntv_district.dat <- list(N=nrow(sntv_district),
                          K=ncol(x.mat.sntv),
                          D=length(unique(sntv_district$distID)),
                          Dist=sntv_district$distID,
                          y=log(sntv_district$invalid),
                          X=x.mat.sntv)

sntv.invalid <- sampling(stanModel.invalid, sntv_district.dat, iter=3000, warmup=1000, 
                         thin=1, chains=3, seed=1234)

						 
						 
				
###
### smdp invalid
###

# election-year matrix
year_list <- sort(unique(smdp_district$year))
year.mat.smdp <- matrix(0, nrow=nrow(smdp_district), ncol=length(year_list)-1)

for(i in 2:length(year_list)){
  year.mat.smdp[smdp_district$year==year_list[i],i-1] <- 1
}

x.mat.smdp <- cbind(smdp_district$kakusu.mean, log(smdp_district$ku_ncand), year.mat.smdp)

smdp_district.dat <- list(N=nrow(smdp_district),
                          K=ncol(x.mat.smdp),
                          D=length(unique(smdp_district$distID)),
                          Dist=smdp_district$distID,
                          y=log(smdp_district$invalid),
                          X=x.mat.smdp)

smdp.invalid <- sampling(stanModel.invalid, smdp_district.dat, iter=3000, warmup=1000, 
                         thin=1, chains=3, 
                         seed=1234)

#save(sntv.invalid, smdp.invalid, file="StanFits_InvalidVote.RData")




###
### summary
###
library(rstan)
#load("StanFits_InvalidVote.RData")

# extract posterior means and 95% credible intervals
invalid.fig <- rbind(summary(sntv.invalid)$summary[2,c(1,4,8)],
                     summary(smdp.invalid)$summary[2,c(1,4,8)])

# Figure 1
#pdf("name_invalid_coef.pdf", height=6, width=6)
par(mar=c(5.1,4.1,2.1,2.1))
plot(1,1,type="n", ylim=c(0.5,2.5), 
     xlim=c(min(invalid.fig[,2])-0.002, max(invalid.fig[,3])+0.002),
     xlab="Posterior Mean and 95% Credible Interval",
     ylab="",
     yaxt="n")
points(invalid.fig[,1], 2:1, pch=20, cex=1.4)
abline(v=0, lty=2, lwd=1)
segments(invalid.fig[,2], 2:1, invalid.fig[,3], 2:1, lwd=2)
segments(invalid.fig[1,2], 2.05, invalid.fig[1,2], 1.95, lwd=2)
segments(invalid.fig[1,3], 2.05, invalid.fig[1,3], 1.95, lwd=2)
segments(invalid.fig[2,2], 1.05, invalid.fig[2,2], 0.95, lwd=2)
segments(invalid.fig[2,3], 1.05, invalid.fig[2,3], 0.95, lwd=2)
axis(2, at=2:1, labels=c("SNTV", "SMDP"), las=1)
#dev.off()


# Table G.2 in the appendix
summary(sntv.invalid, pars=c("beta[1]", "beta[2]", "beta[3]", "sigma_d"))$summary[,c(1,3)]
summary(smdp.invalid, pars=c("beta[1]", "beta[2]", "sigma_d"))$summary[,c(1,3)]